Saline analyzed

install packages

#BiocManager::install("clusterProfiler")
#install.packages("enrichplot")
#installed.packages("patchwork")
#BiocManager::install("SingleR")
#BiocManager::install("celldex")
#BiocManager::install("ensembldb")
#BiocManager::install('limma')
#BiocManager::install("DESeq2")
#organism = "org.Mm.eg.db"
#BiocManager::install(organism)
#BiocManager::install("pathview")
#install.packages("Seurat")
#install.packages("remotes")
#remotes::install_github("satijalab/seurat")
#BiocManager::install("rhdf5")
#install.packages("hdf5r")
#BiocManager::install("NeuCA")
#install.packages('devtools')
#install_github("ggjlab/scMCA")
#install.packages("pheatmap")
#BiocManager::install("scater")
#devtools::install_github("sqjin/CellChat")
#BiocManager::install("ComplexHeatmap")

Load libraries

In this step, we are loading the libraries which are required for analysis and plotting of the data.

library(Seurat)
library(devtools)
library(ggplot2)
library(patchwork)
library(dplyr)
library(clusterProfiler)
library(enrichplot)
# we use ggplot2 to add x axis labels (ex: ridgeplot)
library(stringr)
library(celldex)
library(ensembldb)
library(SingleR)
library(limma)
library(DESeq2)
library(org.Mm.eg.db)
library(rhdf5)
library(hdf5r)
library(pheatmap)
library(ComplexHeatmap)
library(CellChat)
library(patchwork)
options(stringsAsFactors = FALSE)
library(NMF)
library(ggalluvial)

set the colors of the cell types

my_cols <- c('B cells'='#F8766D','DC'='#39568CFF','Endothelial cells'='#CD9600','Epithelial cells'='#00C19A','Fibroblasts'= "#C77CFF",
    'Macrophages'='#00A9FF','Monocytes'='#0CB702','Stromal cells'='#FF61CC', "ILC" = "grey", "Neutrophils" = "Darkred")

load data for mouse 1 and mouse 2

# Mouse 1
saline_mouse_1 <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/Saline_sample_round_one/outs/"

list.files(saline_mouse_1)

##  [1] "analysis"                      "cloupe.cloupe"                
##  [3] "filtered_feature_bc_matrix"    "filtered_feature_bc_matrix.h5"
##  [5] "metrics_summary.csv"           "molecule_info.h5"             
##  [7] "possorted_genome_bam.bam"      "possorted_genome_bam.bam.bai" 
##  [9] "probe_set.csv"                 "raw_feature_bc_matrix"        
## [11] "raw_feature_bc_matrix.h5"      "spatial"                      
## [13] "spatial_enrichment.csv"        "web_summary.html"

saline_mouse_1_data <- Load10X_Spatial(
  saline_mouse_1,
  filename = "filtered_feature_bc_matrix.h5",
  assay = "spatial",
  slice = "saline",
  filter.matrix = TRUE,
  to.upper = FALSE,
  image = NULL)

# Mouse 2
saline_mouse_2 <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/Saline_sample_round_two/outs/"

list.files(saline_mouse_2)

##  [1] "analysis"                      "cloupe.cloupe"                
##  [3] "filtered_feature_bc_matrix"    "filtered_feature_bc_matrix.h5"
##  [5] "metrics_summary.csv"           "molecule_info.h5"             
##  [7] "possorted_genome_bam.bam"      "possorted_genome_bam.bam.bai" 
##  [9] "probe_set.csv"                 "raw_feature_bc_matrix"        
## [11] "raw_feature_bc_matrix.h5"      "saline_web_summary.html"      
## [13] "spatial"                       "spatial_enrichment.csv"

saline_mouse_2_data <- Load10X_Spatial(
 saline_mouse_2,
  filename = "filtered_feature_bc_matrix.h5",
  assay = "spatial",
  slice = "saline",
  filter.matrix = TRUE,
  to.upper = FALSE,
  image = NULL)

data preprocessing

# Mouse 1
saline_mouse_1_data <- SCTransform(saline_mouse_1_data, assay = "spatial", verbose = FALSE)

# Mouse 2
saline_mouse_2_data <- SCTransform(saline_mouse_2_data, assay = "spatial", verbose = FALSE)

Visualization of lymphoid follicle markers

Mouse 1

p1 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Cd19", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p2 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Igha", alpha = c(0.1, 1), pt.size.factor = 1.5)
p3 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Bcl6", alpha = c(0.1, 1), pt.size.factor = 1.5)
p4 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Lta", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p5 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Ltb", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p6 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Cd4", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p7 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Mki67", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p8 <- SpatialFeaturePlot(saline_mouse_1_data, features = "Cxcr5", alpha = c(0.1, 1.5), pt.size.factor = 1.5)

p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8

Mouse 2

p1 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Cd19", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p2 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Igha", alpha = c(0.1, 1), pt.size.factor = 1.5)
p3 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Bcl6", alpha = c(0.1, 1), pt.size.factor = 1.5)
p4 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Lta", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p5 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Ltb", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p6 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Cd4", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p7 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Mki67", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p8 <- SpatialFeaturePlot(saline_mouse_2_data, features = "Cxcr5", alpha = c(0.1, 1.5), pt.size.factor = 1.5)

p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8

dimensionality reduction and clustering

# Mouse 1
saline_mouse_1_analyzed <- RunPCA(saline_mouse_1_data, assay = "SCT", verbose = FALSE)
saline_mouse_1_analyzed <- FindNeighbors(saline_mouse_1_analyzed, reduction = "pca", dims = 1:5)
saline_mouse_1_analyzed <- FindClusters(saline_mouse_1_analyzed, verbose = FALSE, resolution = 0.8)
saline_mouse_1_analyzed <- RunUMAP(saline_mouse_1_analyzed, reduction = "pca", dims = 1:5)

# Mouse 2
saline_mouse_2_analyzed <- RunPCA(saline_mouse_2_data, assay = "SCT", verbose = FALSE)
saline_mouse_2_analyzed <- FindNeighbors(saline_mouse_2_analyzed, reduction = "pca", dims = 1:5)
saline_mouse_2_analyzed <- FindClusters(saline_mouse_2_analyzed, verbose = FALSE, resolution = 0.8)
saline_mouse_2_analyzed <- RunUMAP(saline_mouse_2_analyzed, reduction = "pca", dims = 1:5)

labelling clusters using SingleR pipeline

# First convert the seurat object to single cell experiment object (otherwise the SingleR pipeline will not work)

# Mouse 1
saline_mouse_1_analyzed.sce <- as.SingleCellExperiment(saline_mouse_1_analyzed)

# Mouse 2
saline_mouse_2_analyzed.sce <- as.SingleCellExperiment(saline_mouse_2_analyzed)

#create reference data

ref.data <- celldex::ImmGenData()

# run singleR pipeline to find the cell types

## Mouse 1
cell_types_mouse_1 <- SingleR(test=saline_mouse_1_analyzed.sce, assay.type.test=1, 
    ref=ref.data, labels=ref.data$label.main)

## Mouse 2
cell_types_mouse_2 <- SingleR(test=saline_mouse_2_analyzed.sce, assay.type.test=1, 
    ref=ref.data, labels=ref.data$label.main)

Check if the cell types are labelled with high accuracy

# view cell types

# Mouse 1
table(cell_types_mouse_1$labels)

## 
##           B cells                DC Endothelial cells  Epithelial cells 
##                58                41              1095               262 
##       Fibroblasts               ILC       Macrophages         Monocytes 
##               384                 1              1007                19 
##       Neutrophils     Stromal cells 
##                 2               180

#checking cell types Mouse 1

plotScoreHeatmap(cell_types_mouse_1)

plotDeltaDistribution(cell_types_mouse_1)

summary(is.na(cell_types_mouse_1$pruned.labels))

##    Mode   FALSE 
## logical    3049

# Mouse 2
table(cell_types_mouse_2$labels)

## 
##           B cells                DC Endothelial cells  Epithelial cells 
##                60                11               706                62 
##       Fibroblasts       Macrophages         Monocytes     Stromal cells 
##               162              1534                 9                35

#checking cell types Mouse 2

plotScoreHeatmap(cell_types_mouse_2)

plotDeltaDistribution(cell_types_mouse_2)

summary(is.na(cell_types_mouse_2$pruned.labels))

##    Mode   FALSE    TRUE 
## logical    2570       9

merging metadata of cell types to seurat object

# Mouse 1
saline_mouse_1_analyzed[["ref.data"]] <- cell_types_mouse_1$labels

# Mouse 2
saline_mouse_2_analyzed[["ref.data"]] <- cell_types_mouse_2$labels

Plot umap and spatial plot to localize each cell type

Mouse 1

p1 <- DimPlot(saline_mouse_1_analyzed, group.by = c("ref.data"), reduction = "umap",label = F, pt.size = 1.5, label.size = 5, cols = my_cols)

p2 <- SpatialDimPlot(saline_mouse_1_analyzed, group.by = c("ref.data"), label = FALSE, label.size = 5, cols = my_cols)

p1 + p2

Mouse 2

p1 <- DimPlot(saline_mouse_2_analyzed, group.by = c("ref.data"), reduction = "umap",label = F, pt.size = 1.5, label.size = 5, cols = my_cols)

p2 <- SpatialDimPlot(saline_mouse_2_analyzed, group.by = c("ref.data"), label = FALSE, label.size = 5, cols = my_cols)

p1 + p2

Plot violin plot with top markers associated with B cells and lymphoid follicles

Mouse 1

VlnPlot(saline_mouse_1_analyzed, group.by = c("ref.data"), features = c("Cxcr5", "Mki67", "Ltb", "Cd38", "Ccr6", "Tnfrsf13c"), cols =  my_cols)

Mouse 2

VlnPlot(saline_mouse_2_analyzed, group.by = c("ref.data"), features = c("Cxcr5", "Mki67", "Ltb", "Cd38", "Ccr6", "Tnfrsf13c"), cols =  my_cols)

find markers in each cluster

Mouse 1

B_cell_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data", ident.1 = "B cells", min.pct = 0.25)
head(B_cell_mouse_1, n = 10)

##                  p_val avg_log2FC pct.1 pct.2    p_val_adj
## Ighd      8.255271e-64  1.3507961 0.828 0.118 1.269991e-59
## Cd19      8.357792e-61  1.4494167 0.897 0.155 1.285763e-56
## Il9r      8.629871e-57  0.4171455 0.345 0.018 1.327619e-52
## Ms4a1     2.148405e-52  0.7874505 0.534 0.053 3.305106e-48
## Bank1     3.460125e-52  0.7886827 0.603 0.069 5.323056e-48
## Cr2       3.013761e-50  0.8212669 0.500 0.048 4.636370e-46
## Tnfrsf13c 6.413623e-45  1.0369011 0.707 0.117 9.866717e-41
## Blnk      7.485194e-45  0.9015819 0.741 0.126 1.151522e-40
## Spib      1.319907e-44  1.1623285 0.707 0.118 2.030545e-40
## Fcrla     1.453135e-44  0.6200531 0.483 0.050 2.235503e-40

write.csv(B_cell_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/B_cell_mouse_1.csv")

DC_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "DC", min.pct = 0.25)
head(DC_mouse_1, n = 10)

##               p_val avg_log2FC pct.1 pct.2    p_val_adj
## Zbtb46 2.592652e-08  0.3017862 0.366 0.098 0.0003988535
## Gbp4   3.439860e-08  0.5729455 0.927 0.584 0.0005291880
## Epas1  6.949482e-08 -0.8047154 0.512 0.870 0.0010691084
## Fmo2   1.405798e-07 -0.8504219 0.341 0.683 0.0021626804
## Dffa   1.742615e-07  0.3038923 0.415 0.129 0.0026808388
## Inmt   5.366006e-07 -0.8811450 0.488 0.840 0.0082550640
## Nos2   5.879765e-07  0.5272287 0.902 0.527 0.0090454306
## Rogdi  3.931562e-06  0.3697050 0.634 0.295 0.0604831477
## Lyrm2  4.272549e-06  0.2744930 0.317 0.098 0.0657289010
## Man1c1 4.734111e-06  0.2995503 0.415 0.151 0.0728295623

write.csv(DC_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/DC_mouse_1.csv")

Stromal_cells_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Stromal cells", min.pct = 0.25)
head(Stromal_cells_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Igkc    2.148446e-18 -0.5779135 1.000 1.000 3.305170e-14
## Igha    2.454418e-18 -0.5396479 1.000 0.998 3.775877e-14
## Gsn     2.048482e-15  0.5022576 0.911 0.738 3.151385e-11
## Inmt    8.253679e-14  0.5389418 0.961 0.827 1.269746e-09
## Ly6e    1.938211e-13 -0.3734319 0.956 0.985 2.981744e-09
## Tnfaip2 5.132503e-13 -0.4556705 0.917 0.966 7.895842e-09
## Jchain  1.233062e-12 -0.5806728 0.650 0.864 1.896942e-08
## Ankrd1  3.898149e-10  0.5452627 0.383 0.199 5.996912e-06
## H2-DMb2 1.064259e-09 -0.3793610 0.883 0.931 1.637256e-05
## Ighj1   1.381416e-09 -0.3043841 0.656 0.845 2.125171e-05

write.csv(Stromal_cells_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Stromal_cells_mouse_1.csv")


Endothelial_cells_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Endothelial cells", min.pct = 0.25)
head(Endothelial_cells_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Cldn5   1.027986e-53  0.5591739 0.916 0.772 1.581453e-49
## Igfbp2  3.888415e-48  0.6519764 0.915 0.801 5.981938e-44
## Epas1   8.744963e-45  0.4787279 0.934 0.827 1.345325e-40
## Rasip1  2.717011e-40  0.3987670 0.694 0.464 4.179849e-36
## Tns1    1.906312e-39  0.4303160 0.953 0.872 2.932671e-35
## Timp3   4.363481e-37  0.3786639 0.988 0.942 6.712779e-33
## Clec14a 4.322131e-36  0.3560302 0.554 0.328 6.649166e-32
## Kdr     8.719974e-34  0.3041826 0.526 0.305 1.341481e-29
## Tmcc2   1.355709e-33  0.3613636 0.595 0.375 2.085622e-29
## Igkc    6.254199e-33 -0.3852192 1.000 0.999 9.621460e-29

write.csv(Endothelial_cells_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Endothelial_cells_mouse_1.csv")

Epithelial_cell_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Epithelial cells", min.pct = 0.25)
head(Epithelial_cell_mouse_1, n = 10)

##               p_val avg_log2FC pct.1 pct.2    p_val_adj
## Inmt   1.143840e-19 -0.6655945 0.706 0.847 1.759683e-15
## Cldn5  5.345920e-19 -0.5996306 0.691 0.836 8.224164e-15
## Tns1   4.602813e-16 -0.4876424 0.821 0.908 7.080967e-12
## Ltbp4  5.284907e-16 -0.4639505 0.206 0.480 8.130301e-12
## Timp3  5.455946e-16 -0.4549849 0.905 0.964 8.393427e-12
## Igha   1.112342e-15  0.3166163 1.000 0.998 1.711226e-11
## Igkc   3.078809e-15  0.3298940 1.000 1.000 4.736440e-11
## Pecam1 6.891285e-14 -0.3491782 0.179 0.435 1.060155e-09
## Ramp2  1.604636e-13 -0.3287006 0.092 0.312 2.468573e-09
## Ager   2.040610e-13 -0.3179175 0.973 0.991 3.139274e-09

write.csv(Epithelial_cell_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm+bcg/Epithelial_cell_mouse_1.csv")

Fibroblasts_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Fibroblasts", min.pct = 0.25)
head(Fibroblasts_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Inmt    5.266062e-22  0.4666610 0.966 0.816 8.101309e-18
## Gsn     1.647372e-17  0.4518127 0.878 0.730 2.534317e-13
## Nbl1    1.612737e-14  0.3469147 0.859 0.757 2.481034e-10
## Clec3b  9.523951e-14  0.4607619 0.576 0.420 1.465165e-09
## Fmo2    9.902154e-11  0.2931271 0.794 0.661 1.523347e-06
## Ctss    1.366486e-10 -0.2699602 0.984 0.991 2.102203e-06
## Ltbp4   2.049154e-10  0.3284453 0.581 0.438 3.152419e-06
## Ptprc   1.829641e-09 -0.2610686 0.646 0.768 2.814719e-05
## Rarres2 1.653134e-08  0.3535110 0.729 0.629 2.543181e-04
## Col1a2  1.811325e-08  0.2753930 0.831 0.737 2.786543e-04

  write.csv(Fibroblasts_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Fibroblasts_mouse_1.csv")
  
Macrophages_mouse_1 <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Macrophages", min.pct = 0.25)
  head(Macrophages_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Igha    1.802204e-45  0.3570873 1.000 0.997 2.772510e-41
## Inmt    2.349932e-42 -0.5859870 0.761 0.872 3.615135e-38
## Igfbp2  2.445557e-42 -0.6746716 0.776 0.875 3.762244e-38
## Igkc    3.043211e-37  0.2814937 1.000 1.000 4.681676e-33
## Epas1   4.994383e-35 -0.4661141 0.805 0.895 7.683358e-31
## Cldn5   6.077221e-35 -0.5096620 0.761 0.855 9.349197e-31
## Timp3   6.586934e-33 -0.4033961 0.944 0.966 1.013334e-28
## Tns1    4.589537e-32 -0.4275692 0.865 0.918 7.060544e-28
## Ctss    1.098293e-31  0.3599147 0.998 0.987 1.689614e-27
## Tnfaip2 1.227430e-31  0.3164347 0.995 0.947 1.888278e-27

write.csv(Macrophages_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Macrophages_mouse_1.csv")

Mouse 2

B_cell_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data", ident.1 = "B cells", min.pct = 0.25)
head(B_cell_mouse_2, n = 10)

##                  p_val avg_log2FC pct.1 pct.2    p_val_adj
## Cxcr5     1.397901e-79  1.2684368 0.817 0.089 1.908834e-75
## Bank1     2.629509e-75  1.6224021 0.950 0.146 3.590595e-71
## Cr2       4.492997e-68  1.4266862 0.700 0.077 6.135188e-64
## Il9r      8.464439e-64  0.6025927 0.467 0.031 1.155819e-59
## Pax5      3.205187e-63  1.1707456 0.733 0.090 4.376683e-59
## Fcrl5     5.607384e-63  0.9722059 0.600 0.057 7.656883e-59
## Dtx1      3.669356e-61  0.6956637 0.550 0.047 5.010505e-57
## Ccr6      8.135732e-61  1.9387029 0.950 0.200 1.110934e-56
## Tnfrsf13c 1.292829e-58  1.8747025 0.917 0.178 1.765358e-54
## Ighd      6.070972e-58  1.8733801 0.883 0.167 8.289912e-54

write.csv(B_cell_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/B_cell_mouse_2.csv")

DC_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "DC", min.pct = 0.25)
head(DC_mouse_2, n = 10)

##               p_val avg_log2FC pct.1 pct.2    p_val_adj
## Klhdc9 1.052979e-23  0.3384042 0.273 0.007 1.437843e-19
## Clec9a 5.968194e-19  0.4251603 0.364 0.016 8.149569e-15
## Cfap74 8.906069e-19  0.5914211 0.455 0.025 1.216124e-14
## Mst1   1.538070e-14  0.3311671 0.273 0.012 2.100234e-10
## Zfp160 1.356914e-13  0.4892451 0.455 0.035 1.852866e-09
## Erich3 1.016938e-11  1.0236002 0.545 0.059 1.388628e-07
## Lrrc43 4.023684e-11  0.9656281 0.636 0.086 5.494340e-07
## March4 5.542413e-11  0.4979459 0.364 0.028 7.568166e-07
## Fzd8   1.651961e-10  0.3223094 0.273 0.017 2.255753e-06
## Gp2    6.170181e-10  0.4113967 0.273 0.018 8.425383e-06

write.csv(DC_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/DC_mouse_2.csv")

Stromal_cells_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Stromal cells", min.pct = 0.25)
head(Stromal_cells_mouse_2, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Myh11   2.799853e-12  0.9352038 0.571 0.163 3.823199e-08
## Cst3    4.295150e-12  0.6681956 1.000 1.000 5.865027e-08
## Des     7.746254e-12  1.0718122 0.857 0.488 1.057751e-07
## Tagln   7.373254e-11  0.9235898 0.771 0.317 1.006818e-06
## Nbl1    4.602124e-10  0.7027256 1.000 0.985 6.284200e-06
## Nphp4   5.992370e-10  0.3655504 0.314 0.060 8.182581e-06
## Gsn     8.847970e-10  0.8142324 1.000 0.927 1.208190e-05
## Fth1    8.466740e-09 -0.4846324 1.000 1.000 1.156133e-04
## Pcp4l1  9.268633e-09  0.4812146 0.457 0.126 1.265632e-04
## Col14a1 1.227715e-08  0.3198091 0.314 0.066 1.676445e-04

write.csv(Stromal_cells_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Stromal_cells_mouse_2.csv")


Endothelial_cells_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Endothelial cells", min.pct = 0.25)
head(Endothelial_cells_mouse_2, n = 10)

##               p_val avg_log2FC pct.1 pct.2    p_val_adj
## Cldn5  1.411386e-98  0.8332450 0.986 0.896 1.927248e-94
## Timp3  2.632595e-76  0.5422023 1.000 0.996 3.594809e-72
## Epas1  1.200408e-75  0.6283701 0.992 0.953 1.639157e-71
## Acer2  2.821727e-67  0.5590838 0.993 0.952 3.853068e-63
## Igfbp2 4.746039e-67  1.0716958 0.877 0.690 6.480717e-63
## Egfl7  1.783598e-66  0.5758111 0.992 0.915 2.435503e-62
## Ace    4.473554e-62  0.6363960 0.967 0.874 6.108638e-58
## Tspan7 6.054914e-62  0.6679605 0.856 0.612 8.267984e-58
## Plvap  1.008573e-58  0.4802857 0.997 0.981 1.377207e-54
## Pecam1 3.042090e-56  0.5633490 0.939 0.798 4.153974e-52

write.csv(Endothelial_cells_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Endothelial_cells_mouse_2.csv")

Epithelial_cell_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Epithelial cells", min.pct = 0.25)
head(Epithelial_cell_mouse_2, n = 10)

##                      p_val avg_log2FC pct.1 pct.2    p_val_adj
## 1700016K19Rik 4.489550e-59  0.9366679 0.629 0.068 6.130480e-55
## Dynlrb2       5.376514e-48  0.7171786 0.484 0.048 7.341630e-44
## 1700013F07Rik 1.266140e-47  0.5067919 0.387 0.030 1.728915e-43
## Wdr63         1.681633e-45  0.6232098 0.484 0.050 2.296271e-41
## Kndc1         4.331139e-44  0.6357987 0.468 0.049 5.914171e-40
## March4        7.567509e-43  0.4041338 0.323 0.023 1.033343e-38
## Tekt1         2.856782e-41  1.0867857 0.629 0.101 3.900935e-37
## Rp1           3.613835e-41  1.0314188 0.629 0.100 4.934692e-37
## Alox12e       5.386355e-41  0.3617788 0.274 0.017 7.355067e-37
## Dnah5         1.044865e-39  0.7036400 0.532 0.072 1.426763e-35

write.csv(Epithelial_cell_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Epithelial_cell_mouse_2.csv")

Fibroblasts_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Fibroblasts", min.pct = 0.25)
head(Fibroblasts_mouse_2, n = 10)

##               p_val avg_log2FC pct.1 pct.2    p_val_adj
## Col1a2 3.272175e-17  0.4468853 1.000 0.998 4.468154e-13
## Igha   5.584053e-15  1.5201650 0.994 1.000 7.625025e-11
## C3     1.046965e-14  0.2901238 1.000 1.000 1.429631e-10
## Apoe   6.533885e-14  0.4792034 0.957 0.892 8.922020e-10
## Gbp2b  2.702983e-13 -0.3358854 1.000 1.000 3.690923e-09
## Igkc   6.114446e-13  1.3331430 0.994 1.000 8.349276e-09
## Col1a1 6.408404e-13  0.4230387 0.988 0.969 8.750676e-09
## Ctss   9.100994e-13 -0.3472002 1.000 1.000 1.242741e-08
## Mzb1   2.642364e-12  0.7704413 0.691 0.467 3.608148e-08
## Col3a1 5.232022e-12  0.4323003 1.000 0.994 7.144325e-08

  write.csv(Fibroblasts_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Fibroblasts_mouse_2.csv")
  
Macrophages_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Macrophages", min.pct = 0.25)
  head(Macrophages_mouse_2, n = 10)

##               p_val avg_log2FC pct.1 pct.2    p_val_adj
## Ctss   9.097211e-77  0.4498103 1.000 1.000 1.242224e-72
## Fth1   7.823133e-71  0.3583520 1.000 1.000 1.068249e-66
## Psap   5.709096e-56  0.3176438 1.000 1.000 7.795771e-52
## Cfb    1.370117e-55  0.4771078 0.994 0.988 1.870895e-51
## Itgb2  4.442422e-51  0.4104689 0.998 0.984 6.066127e-47
## Il1rn  1.955168e-49  0.5930291 0.984 0.958 2.669781e-45
## Lgals3 3.729628e-48  0.4458466 0.999 0.995 5.092807e-44
## Nos2   1.335862e-47  0.5751555 0.977 0.939 1.824119e-43
## Cyba   2.697403e-46  0.2925800 1.000 1.000 3.683303e-42
## Ctsb   7.022225e-42  0.3247272 1.000 1.000 9.588848e-38

write.csv(Macrophages_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Macrophages_mouse_2.csv")

Monocytes_mouse_2 <- FindMarkers(saline_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Monocytes", min.pct = 0.25)
  head(Monocytes_mouse_2, n = 10)

##               p_val avg_log2FC pct.1 pct.2    p_val_adj
## Nog    9.502314e-11  0.3833897 0.333 0.020 1.297541e-06
## Nuf2   3.420788e-08  0.3746294 0.333 0.027 4.671086e-04
## Espl1  3.704791e-07  0.3686375 0.333 0.032 5.058893e-03
## Cd5    1.054815e-06  0.3659221 0.333 0.034 1.440349e-02
## Sfmbt2 1.310108e-06  0.3626703 0.333 0.034 1.788952e-02
## Crmp1  1.440750e-06  0.3653796 0.333 0.035 1.967344e-02
## Wbp1   2.018414e-06  0.4423256 0.444 0.061 2.756144e-02
## Zfp930 8.394765e-06  0.4700501 0.333 0.040 1.146305e-01
## Cbr4   2.050164e-05  0.4826788 0.556 0.108 2.799499e-01
## Camk4  5.347287e-05  0.4656349 0.556 0.116 7.301720e-01

write.csv(Monocytes_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Monocytes_mouse_2.csv")

Macrophage_cell_deseq <- FindMarkers(saline_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Macrophages", test.use = "DESeq2")
write.csv(Macrophage_cell_deseq, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Macrophage_cell_deseq.csv")

Cell-cell interaction using cell chat

saline_cellchat <- createCellChat(object = saline_mouse_1_analyzed, group.by = "ref.data")

## The `data` slot in the default assay is used. The default assay is SCT 
## The `meta.data` slot in the Seurat object is used as cell meta information 
## The cell groups used for CellChat analysis are  B cells DC Endothelial cells Epithelial cells Fibroblasts ILC Macrophages Monocytes Neutrophils Stromal cells

CellChatDB <- CellChatDB.mouse
showDatabaseCategory(CellChatDB)

CellChatDB.use <- CellChatDB
saline_cellchat@DB <- CellChatDB.use


saline_cellchat <- subsetData(saline_cellchat)
saline_cellchat <- identifyOverExpressedGenes(saline_cellchat)
saline_cellchat <- identifyOverExpressedInteractions(saline_cellchat)


saline_cellchat <- computeCommunProb(saline_cellchat)
saline_cellchat <- filterCommunication(saline_cellchat, min.cells = 10)

## The cell-cell communication related with the following cell groups are excluded due to the few number of cells:  ILC Neutrophils

saline_cellchat <- computeCommunProbPathway(saline_cellchat)
saline_cellchat <- aggregateNet(saline_cellchat)

groupSize <- as.numeric(table(saline_cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(saline_cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions")
netVisual_circle(saline_cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")

# See interaction of each cell type with the another

mat1 <- saline_cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat1)) {
  mat2 <- matrix(0, nrow = nrow(mat1), ncol = ncol(mat1), dimnames = dimnames(mat1))
  mat2[i, ] <- mat1[i, ]
  netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat1), title.name = rownames(mat1)[i])
}

##Visualize signalling pathways in each cell

# Identify signaling pathways showing significant communications 
saline_cellchat@netP$pathways

##  [1] "COLLAGEN"    "APP"         "MIF"         "THBS"        "VEGF"       
##  [6] "FN1"         "LAMININ"     "MHC-II"      "ICAM"        "ITGAL-ITGB2"
## [11] "COMPLEMENT"  "PECAM1"      "CD52"        "CD22"        "CD45"       
## [16] "JAM"         "PARs"        "CHEMERIN"    "SEMA4"       "CXCL"       
## [21] "EPHA"        "PDGF"        "GAS"         "HSPG"        "FGF"        
## [26] "SEMA3"       "AGRN"        "PTPRM"       "CDH5"        "SEMA7"      
## [31] "NOTCH"       "SEMA6"       "IL16"

# Evaluate one of the pathways
pathways.show <- c("NOTCH") 
vertex.receiver = seq(2,5) # a numeric vector. 
?netVisual_aggregate(saline_cellchat, signaling = pathways.show,  vertex.receiver = vertex.receiver)

# Chord diagram
par(mfrow=c(1,1))
netVisual_aggregate(saline_cellchat, signaling = pathways.show, layout = "chord")

# Heatmap
par(mfrow=c(1,1))
netVisual_heatmap(saline_cellchat, signaling = pathways.show, color.heatmap = "Reds")

Compute the contribution of each ligand-receptor pair to the overall signaling pathway

netAnalysis_contribution(saline_cellchat, signaling = pathways.show)

pairLR.MHC <- extractEnrichedLR(saline_cellchat, signaling = pathways.show, geneLR.return = FALSE)
LR.show <- pairLR.MHC[1,] # show one ligand-receptor pair

# Hierarchy plot
vertex.receiver = seq(1,4) # a numeric vector
netVisual_individual(saline_cellchat, signaling = pathways.show,  pairLR.use = LR.show, vertex.receiver = vertex.receiver)

## [[1]]

Visualize cell-cell communication mediated by multiple signaling pathways

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_bubble(saline_cellchat, sources.use = 1, targets.use = c(2:8), remove.isolate = FALSE)

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# show all the interactions sending from Inflam.FIB
netVisual_chord_gene(saline_cellchat, sources.use = 1, targets.use = c(2:8), lab.cex = 0.5,legend.pos.y = 30)

Compute the network centrality scores

saline_cellchat <- netAnalysis_computeCentrality(saline_cellchat, slot.name = "netP") 

# the slot 'netP' means the inferred intercellular communication network of signaling pathways

# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
netAnalysis_signalingRole_network(saline_cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)

Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways

ht1 <- netAnalysis_signalingRole_heatmap(saline_cellchat, pattern = "outgoing")
ht2 <- netAnalysis_signalingRole_heatmap(saline_cellchat, pattern = "incoming")
ht1 + ht2

Identify and visualize communication patterns (outgoing signal)

selectK(saline_cellchat, pattern = "outgoing")

nPatterns = 3
saline_cellchat <- identifyCommunicationPatterns(saline_cellchat, pattern = "outgoing", k = nPatterns)

# river plot
netAnalysis_river(saline_cellchat, pattern = "outgoing")

# dot plot
netAnalysis_dot(saline_cellchat, pattern = "outgoing")

Identify and visualize communication patterns (incoming signal)

selectK(saline_cellchat, pattern = "incoming")

nPatterns = 3
saline_cellchat <- identifyCommunicationPatterns(saline_cellchat, pattern = "incoming", k = nPatterns)

# river plot
netAnalysis_river(saline_cellchat, pattern = "incoming")

# dot plot
netAnalysis_dot(saline_cellchat, pattern = "incoming")

saline_cellchat <- computeNetSimilarity(saline_cellchat, type = "functional")
saline_cellchat <- netEmbedding(saline_cellchat, type = "functional")

## Manifold learning of the signaling networks for a single dataset

#> Manifold learning of the signaling networks for a single dataset
saline_cellchat <- netClustering(saline_cellchat, type = "functional")

## Classification learning of the signaling networks for a single dataset

#> Classification learning of the signaling networks for a single dataset
# Visualization in 2D-space
netVisual_embedding(saline_cellchat, type = "functional", label.size = 3.5)

Gene set enrichment

organism = org.Mm.eg.db

prepare input data for gene set enrichment

# reading in data from deseq2
df = read.csv("/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/saline/Macrophage_cell_deseq.csv", header=TRUE)

# we want the log2 fold change 
original_gene_list <- df$avg_log2FC

# name the vector
names(original_gene_list) <- df$X

# omit any NA values 
gene_list<-na.omit(original_gene_list)

# sort the list in decreasing order (required for clusterProfiler)
gene_list = sort(gene_list, decreasing = TRUE)

Gene set enrichment

#gene set enrichment function for B cell cluster

gse <- gseGO(geneList=gene_list, 
             ont ="ALL", 
             keyType = "SYMBOL", 
             minGSSize = 3, 
             maxGSSize = 800, 
             pvalueCutoff = 0.05, 
             verbose = TRUE, 
             OrgDb = organism, 
             pAdjustMethod = "none",
             eps = 0.0)

plot graph for upregulated and downregulated genes

require(DOSE)
gse_dataframe <- as.data.frame(gse)

## count the gene number
gene_count<- gse@result %>% group_by(ID) %>% summarise(count = sum(str_count(core_enrichment, "/")) + 1)

## merge with the original dataframe
gse_dataframe<- left_join(gse@result, gene_count, by = "ID") %>% mutate(GeneRatio = count/setSize)

write.csv(gse_dataframe, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm+bcg/gse_dataframe.csv")

ggplot(gse_dataframe[c(3,9,17,21,27,29,41,58,59,63,89,102,109,118,136,139,145,163),],
             aes(x = GeneRatio, y = Description)) + 
             geom_point(aes(size = GeneRatio, color = p.adjust)) +
             theme_bw(base_size = 14) +
             scale_colour_gradient(limits=c(0, 0.10), low="red") +
             ylab(NULL) +
             ggtitle("GO pathway enrichment")